# Required packages
require(tidyverse)
require(sf)
require(spdep)
require(spatialreg)
require(lme4)
require(lmerTest)

# Version used:
# sf_1.0-9
# spdep_1.2-5
# spatialreg_1.2-5 
# lme4_1.1-30
# lme4_1.1-30

# Load data
load("01_Data/spatial_dimension_data.Rdata")

#### Models ####

#### 1 - Linear models #### 

# NOTE ON VARIABLES
# delay = deviation from legal due date
# time = drive time in hours
# risk_factor = risk level

# Pooled model on location level
lm0 <- lm(delay ~ time + risk, location_level)

# Models on plant level
lm1 <- lm(delay ~ time + 
            risk_factor, 
            model_base)

lm2 <- lm(delay ~
            time +
            risk_factor +
            pressure +
            agent +
            agencies +
            sector_1+
            sector_2+
            sector_3+
            year,
          model_base)

lm3 <- lm(interval ~
            time +
            risk_factor +
            pressure +
            agent +
            agencies +
            sector_1+
            sector_2+
            sector_3+
            year,
          model_base)

#### 2 - Calculate weight matrix ####

# Look for critical threshold
k1 <- knn2nb(knearneigh(location_coordinates, k = 1, longlat = T))
critical.threshold <- max(unlist(nbdists(k1, location_coordinates, longlat = T)))

# Set distance band to criitcal threshold
nb.dist.band <-
  dnearneigh(location_coordinates, 0, critical.threshold, longlat = T)
distances <- nbdists(nb.dist.band, location_coordinates, longlat = T)
invd1 <- lapply(distances, function(x)
  (1 / x))

invd.weights <- nb2listw(nb.dist.band, glist = invd1, style = "W")

# Expand matrix to fit full data set
mat2 <- as.data.frame(spdep::listw2mat(invd.weights))
mat2 <- as.data.frame(lapply(mat2, rep, location_level$n))
mat2 <- as.data.frame(t(mat2))
mat2 <- as.data.frame(lapply(mat2, rep, location_level$n))
mat2 <- t(mat2)


invd.weights_large <- spdep::mat2listw(mat2)
rm(distances, invd1, k1, mat2, nb.dist.band)

# save weight matrix
save(invd.weights,
     invd.weights_large,
     file = "01_Data/weight_matrix.Rdata")

### Test for resiudal spatial autocorrelation ###
lm.morantest(lm0, listw = invd.weights)
lm.morantest(lm1, listw = invd.weights_large)
lm.morantest(lm2, listw = invd.weights_large)
lm.morantest(lm3, listw = invd.weights_large)

#### 3 - Spatial autoregressive models ####
sm1 <- lagsarlm(delay ~
                  time +
                  risk_factor,
                data = model_base,
                listw = invd.weights_large)

impacts(sm1, listw = invd.weights_large)

sm2 <- spatialreg::lagsarlm(
  delay ~ time +
    risk_factor +
    pressure +
    agent +
    agencies +
    sector_1+
    sector_2+
    sector_3+
    year,
  data = model_base,
  listw = invd.weights_large
)

impacts(sm2, listw = invd.weights_large)

sm3 <- lagsarlm(
  interval ~ time +
    risk_factor +
    pressure+
    agent +
    agencies +
    sector_1+
    sector_2+
    sector_3+
    year,
  data = model_base,
  listw = invd.weights_large
)

impacts(sm3, listw = invd.weights_large)

# Spatial autoregressive models, pooled
sm4 <- lagsarlm(delay ~ time + risk,
                data = location_level,
                listw = invd.weights)

impacts(sm4, listw = invd.weights)



#### 4 - Mixed effect models #### 

model_base$district[model_base$district == 1] <- "Stuttgart"
model_base$district[model_base$district == 2] <- "Karlsruhe"
model_base$district[model_base$district == 3] <- "Freiburg"
model_base$district[model_base$district == 4] <- "Tübingen"

lmer1 <- lmer( delay ~ time +
                 risk_factor +
                 (1|qtr:district)+
                 (1|id),
               REML = F,
               model_base)


lmer2 <- lmer( delay ~ time +
                 risk_factor +
                 pressure +
                 agent +
                 agencies +
                 sector_1 +
                 sector_2 +
                 sector_3 +
                 (1|id)+
                 (1|qtr:district),
               REML=F,
               model_base)

lmer3 <- lmer( interval ~ time +
                 risk_factor +
                 pressure +
                 agent +
                 agencies +
                 sector_1 +
                 sector_2 +
                 sector_3 +
                 (1|id)+
                 (1|qtr:district),
               REML=F,
               model_base)

#### 5 - Save all models ####
save(lm1,
     lm2,
     lm3,
     sm1,
     sm2,
     sm3,
     lmer1,
     lmer2,
     lmer3,
     file="01_Data/spatial_dimension_models.Rdata"
)

save(invd.weights,
     invd.weights_large,
     file="01_Data/spatial_dimension_weight_matrix.Rdata"
)

rm(list = ls())


